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A robust and efficient frequency dependent and non-local exchange-correlation f xc (r, r'; o») is 
derived by imposing time-dependent density-functional theory (TDDFT) to reproduce the many- 
body diagrammatic expansion of the Bethe-Salpeter polarization function. As an illustration, we 
compute the optical spectra of LiF, S1O2 and diamond and the finite momentum transfer energy-loss 
spectrum of LiF. The TDDFT results reproduce extremely well the excitonic effects embodied in 
the Bethe-Salpeter approach, both for strongly bound and resonant excitons. We provide a working 
expression for f xc that is fast to evaluate and easy to implement. 
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Since the 30's, excitons have been ubiquitous in our 
understanding of the optics of bulk materials, surfaces, 
nanostructures and organic/bio-molecules 0. Only re- 
cently, however, has the first principle description of ex- 
citons in the optical absorption of extended systems been 
achieved, by solving the Bethe-Salpeter equation (BSE) 
of Many-Body Perturbation Theory (MBPT) Q. The 
solution of the BSE is usually cast into an equivalent 
Hamiltonian problem whose dimension increases with the 
number of k-points and number of valence and conduc- 
tion bands. However, even if the BSE results reproduce 
well the experimental spectra for semiconductors and in- 
sulators, the complexity of the calculations impedes the 
application of this technique to large systems such as 
nanostructures and complex surfaces. 

An alternative approach to the study of correlation in 
many-body systems is given by density-functional the- 
ory, in its static (DFT) [3| and time dependent formula- 
tions (TDDFT) Q. Similar to the paradigm of DFT for 
ground-state properties, TDDFT has emerged as a very 
powerful tool for the description of excited states. In 
principle TDDFT is exact for neutral excited-state prop- 
erties, and its simplicity relies on the fact that two-point 
response functions are needed instead of the four-point 
function of the BSE 1]. TDDFT casts all many-body 
effects into the dynamical exchange-correlation kernel 
f xc (r,t,r';t') = 5v xc (r, t) /5p (r', t'), where v xc (r,t) is 
the corresponding time-dependent exchange-correlation 
potential. It was early recognized || that, in extended 
systems, the standard approximations for v xc - local den- 
sity (LDA) or generalized gradient (GGA) - that work ex- 
tremely well for ground state properties, fail to describe, 
among other effects, the band-gap of insulators and semi- 
conductors and the excitonic effects in the optical and 
energy-loss spectra |l| . Recently, promising results have 
been obtained within a polarization dependent functional 
derived in the framework of current-DFT 6] and within 



the exact-exchange DFT approach Q. The calculated 
optical spectrum of silicon exhibits excitonic effects in 
qualitative agreement with experiment. However, em- 
pirical cutoffs were introduced to construct f xc 0, that 
somehow account for the screening of the electron-hole 
interaction. Furthermore, to-date the calculations of the 
absorption spectra of solids beyond time-dependent LDA 
were performed in semiconductors characterized by weak 
continuum excitonic effects 0, H IS EU EH 113 • Conse- 
quently it remains open whether or not strong electron- 
hole effects (e.g. bound excitons) in the optical and 
energy-loss spectra can be described within TDDFT. 
The present letter resolves unambiguously this contro- 
versy, by deriving a novel f xc that indeed accounts for 
excitonic effects in semiconductors and wide-gap insula- 
tors. 

The challenge is how to construct f xc . To address this 
goal we benefit from the good performance of MBPT 
response functions and build an f xc that mimics those 
results. This f xc is derived by imposing TDDFT to re- 
produce the perturbative expansion of the BSE in terms 
of the screened Coulomb interaction at any order. Our 
derivation starts from the usual TDDFT equation for the 
irreducible response function P (q, uj) Q 

P (q, w) = P (0) (q, uj) + P (0) (q, w) f xc (q, w) P (q, w) , 

(1) 

where P^ (q, oj) is the "exact" Kohn-Sham DFT re- 
sponse function for momentum transfer q. All quantities 
are two-point functions (matrices in reciprocal space). 
The microscopic dielectric matrix is obtained from P and 
the Coulomb potential v as e (q, to) = 1 vP (q, uj). Now 
we make the connection with MBPT. First, P^ (q, uj) is 
approximated by the independent-quasiparticle (QP) re- 
sponse, calculated in the GW scheme [bj. Second, we 
assume that there exists an f xc (q, uS) that reproduces the 
BSE spectra i.e., we impose P (q,cj) = Pbse (q,w). 
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Without loss of generality, we restrict our derivation to 
the resonant part of P (q, to): 



Pi 



Gi,G. 



(q,w) 



2i 

n 



Y, ^ 1 (q,Gi)5 Klj K 2 (q,a;)$K 2 (q,G 2 ), (2) 



K,.K 2 



where f2 the crystal volume, and K = (c v k) a gen- 
eralized index to describe the space of electron-hole 
states. The oscillators <&k are given by $k (q, G) = 
(ck|e^ q+G - ) ' r |t;k — q), in terms of the conduction and va- 
lence Kohn-Sham states. In Eq. J2J S (q, u>) is the solu- 
tion of the BSE 0] 

S (q, w) = SS^ (q, w ) + SS(0) ( q , w ) W (q) S (q, w) , 



We illustrate the reliability of the present TDDFT ap- 
proach in three prototype systems: LiF, Si0 2 and dia- 
mond 01 • I n these three systems the role of excitonic ef- 
fects in the optical spectrum and EELS has been already 
analyzed within the BSE 0, ^| . Si0 2 is characterized 
by four strong excitonic peaks at 10.3, 11.3, 13.5 and 17.5 
eV, none of them below the QP gap of 10.1 eV, except 
for a bound triplet exciton optically inactive. Moreover, 
the exciton at 10.3 eV corresponds to a strongly corre- 
lated resonant state with a large degree of spatial local- 
ization (2-3 bond lengths) 0] . The spectrum of LiF is 
dominated by a strongly bound exciton (~ 3 eV binding 
energy) jl6j . Last, in diamond, the electron-hole interac- 
tion produces a drastic modification of the independent 
QP spectrum by shifting optical oscillator strength from 
high to low energies. 
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in terms of the electron-hole 



QP energies (e^ p ). W (q) is the Coulombic part of 
the Bethe-Salpeter kernel, WkjKj (q) = *(ciki, i^ka — 
q\W (ri, r 2 ) |c 2 k 2 , Viki — q), with W(ri,r 2 ) the stati- 
cally screened electron-hole interaction. Eq. (JIJ can be 
transformed into an equation for f xc (q, u>) by putting 
P(q,o;)=p( )(q,w) + 5P(q,cj) 

P (0) (q^)fxc (q,^)P (0) (q,co) = 

6P (q, to) - SP (q, to) f xc (q, «) P (0) (q, to) . (4) 

The advantage of Eq. (0J is that, by expanding SP (q, uf) 

- (n) 

in terms of W, SP (q, w) = 5P (q, w), it is possi- 
ble to write the n-th order (fie"* (q, w)) contribution to 
fee (q, ui) in an iterative form: 
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with SP (q,w) = 0. This ends our derivation. 

Other attempts to build from MBPT rely on ei- 
ther fully solving the BSE [lflj or by imposing f xc to be 
static and linear in W 0. However, the static constraint 
has to be released in the practical solution for the optical 
spectra of silicon and SiC [l^. This internal inconsis- 
tency is solved in the present work as Eq. JSJ) provides 
a systematic and consistent treatment of the frequency 
dependence and local field effects (LFE) of f xc (q, u>) at 
the cost of making f xc (q, uS) a non-linear functional of 

W . The scheme proposed in Refs. guana appears 

naturally as an approximate solution of Eq. JSJ . 
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FIG. 1: Optical absorption spectra calculated within the 
BSE (black line), an d T DDFT (red line) compared with exper- 
iments (blue circles) |l7l and with the independent-QP calcula- 
tion (dashed line). The results obtained using a scalar f xc (green 
line) are shown to stress the importance of the LFE in f xc . The fre- 
quency dependence of the head of t xc , lim q _,o K [g 2 fuc (q, iif)] 
is shown in the insets (note the different scale for Si02 , LiF and 
C). 



The TDDFT results for the optical absorption of those 
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three elements using the first order f xc of Eq. (0 (f : 
[p(0)] -1 £P (1) [p(°)] _1 ) are compared in Figured- 
line) with the BSE spectra (black line) and experi- 
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ments (blue circles) . To put in evidence the strong 
role played by excitonic effects, the independent-QP 
spectrum is also shown (dashed line). In all three cases 
the agreement between TDDFT and BSE is excellent. 
This good performance of holds not only for the head 
of e (q, ui) but also for the off-diagonal elements. 

The computational cost of our scheme is mainly given 
by the size of e (q, oj) that is proportional to the degree 
of inhomogeneity of the induced density, i.e. it is re- 
lated to the microscopic LFE. In the present TDDFT 
calculations of the response functions we have two con- 
tributions to LFE: the standard contribution from the 
Hartree potential as included in a calculation within the 
random-phase approximation (RPA) and the one com- 
ing from the exchange-correlation effects embodied in f xc 
(e.g., electron-hole attraction). In all the systems we have 
looked at, it turns out that the size of the LFE from f xc 
is smaller or equal to that of the Hartree contribution. 
This general behavior of f xc in extended systems is due 
to the spatial localization of the excitonic wavefunctions, 
that exceeds the typical length scale for the density varia- 
tions (bond-length) . As the proper description of LFE is 
achieved with matrix sizes much smaller than the typical 



size of the BSE Hamiltonian (of the order of 10 
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the present TDDFT approach should be computationally 
favorable with respect to BSE. 

To elucidate the role of the LFE in f xc we show in Fig.Q] 
with green-line the results imposing a first -order scalar, 

but frequency dependent f xc : f£) (q — > 0, ui) 

L Jg=g'=o 
(whose real part is shown in the insets of Fig.^). For di- 
amond, as well as other semiconductors not shown here, 
this scalar f x x J works very well (see Fig.^). In particular, 
the head of f xc is a smooth function of frequency below 
the QP-gap. This result supports the long-range model 
of Ref. 9]. However, SiC>2 and LiF are not at all well de- 
scribed by this scalar f xc . Only the main peak is qual- 
itatively reproduced but not at the correct energy, and 
even unphysical regions of negative absorption appear for 
LiF at high energies. A good spectrum is only achieved 
when the f x V matrix dimension is set to 267 x 267 and 
59 x 59 for SiC>2 and LiF, respectively. Increasing the f xc 
matrix size does not introduce any change in the spec- 
trum. Furthermore, the head of f xc is strongly frequency 
dependent in order to describe the high-energy features 
of the spectra (see inset in Fig.^). In conclusion, for the 
description of the absorption and energy-loss spectra in 
systems with continuum excitonic effects (e.g. diamond) 
the frequency dependence and the microscopic LFE of f xc 
are not important while they become crucial in wide-gap 
insulators with bound- excitons (e.g. LiF) or strongly cor- 
related resonant states (e.g. SiC>2 )■ 

To show the robustness and transferability of the per- 
turbative approach to f xc we calculated the EELS of LiF 
for a finite transfer momentum q along the TX direc- 
tion where previous BSE calculations and experimental 



results are available 0]. This is a stringent test as the 
description of EELS needs causal response functions, in- 
cluding the anti-resonant part 0] . In the present theory 
a causal f xc can be easily obtained by inserting in Eq. Q 
the causal P (q, uo) obtained from the time-ordered P 
of MBPT. The results of this calculation are presented 
in Fig-EI f° r a I s * order (dot-dashed line) and a 2 nd or- 
der (dashed line) causal f xc . In contrast to the results 
of Fig.Q] for the optical absorption, we need to go to a 
second order causal i xc to restore the good agreement 
between TDDFT and BSE. Still, a first order (q,w) 
gives very reasonable EEL spectrum. The results for the 
optical absorption and EEL highlight that the first order 
fxc (q, w) embodies all relevant many body effects at the 
BSE level even though it is a contracted two-point func- 
tion. Moreover higher order contributions to f xc yield 
minor modifications to the spectrum because strong can- 
cellations occur at any order, except for the first, of the 
perturbative expansion of f xc (q, lu) [Eq. JSJ]. 
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FIG. 2: (Color online) Calculated EELS of LiF for a momentum 
transfer of q = 0.51" X: BSE (full line), 1 st order i xc (dot-dashed 
line), 2 nd order f xc (dashed line), and ind ependent-QP (dots). Ex- 
periment (circles) are taken from Ref. llSl . 



A numerical remark about the results in Figs.n[21 is 
relevant now. The direct application of Eq. (J5J) for f xc 
leads to spurious oscillations in the calculated optical 
spectra. Those oscillations are moderate in the case of di- 
amond but they tend to destroy the spectra for the case of 
LiF and SiC>2 ■ The intensity of the oscillations increases 
with the order of fi"'* and eventually, gives rise to non- 
physical regions of negative absorption. The reason for 
this numerical pathology stems from the way the Bethe- 
Salpeter kernel acts on the spectra: (i) redistributing the 
optical oscillator strength, (ii) shifting rigidly the spectra 
to account for the diagonal of the Bethe-Salpeter kernel 
A q = Wk,k (q), that should vanish in the limit of infinite 
k-point sampling. However with a finite k-point grid, 
corresponding to a fully converged BSE spectrum, wc 



get A c 



-0.3 eV in diamond, and A c 



-0.9 eV 



in LiF and SiC>2 ■ The diagonal of W (q) appears in 
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f xc (q, uj) through a series expansion in A q |20j that is 
meaningful only when A q is sufficiently small. As this 
is not the case for SiC>2 and LiF, natural oscillations are 
found in the naive application of the perturbative expan- 
sion of Eq. @ . To circumvent this issue we included the 
diagonal part of W (q) into the independent QP response 

function p(°) (q,w) and let fj c (Qj^) account explicitly 
for the off-diagonal contributions to the Bethe-Salpeter 
kernel. Using this idea the higher order corrections to 
%xc (q, are not only well defined but numerically sta- 
ble at all orders with the same k-point sampling used in 
a standard BSE calculation. 

In conclusion, we have shown that TDDFT with a 
first order fic (q, u>) reproduces the optical and energy- 
loss spectra for a large class of materials: insulators and 
wide-gap insulators. In particular, bound excitons are 
described within TDDFT. Still, the direct implementa- 
tion of fjc (q, to) is cumbersome. However, by looking 

at the analytic properties of 5P^ ^ (q, u>) (for an explicit 
expression see Ref. 01) we can single out the contribu- 
tion of degenerate non-interacting electron-hole states in 
Eq. J2J and write a general expression for fJe (q, ui): 



fW (q, W ) = £[p<°J (q,u/) 
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(6) 



Here u>' — 
independent 



A q and the sum runs to all 



electron-hole 
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degenerate states. Eq. © is the main practical result 
of this letter. It is very fast to compute [2l| as it has 
the form of a non-interacting polarization function with 
modified residuals (Q, R) that are evaluated only once 
as a result of two simple matrix-vector multiplications. 
Also, Eq. © can be made causal and be extended to 
higher orders of the perturbative expansion of f xc (q, u>). 
This implementation opens the way for calculations 
of the response function of nanostructures and low- 
dimensional systems within TDDFT. Work along this 
line is under progress. 
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holds for the description of QP-shifts. the direct diagonalization of the BSE hamiltonian. 

In the case of diamond the solution of the TDDFT equa- 
tion using Eq. © is two orders of magnitude faster than 



